Determination of the diffusion constant using phase-sensitive measurements 
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We apply a pulsed-light interferometer to measure both the intensity and the phase of light that is 
transmitted through a strongly scattering disordered material. From a single set of measurements we 
obtain the time-resolved intensity, frequency correlations and statistical phase information simulta- 
neously. We compare several independent techniques of measuring the diffusion constant for diffuse 
propagation of light. By comparing these independent measurements, we obtain experimental proof 
of the consistency of the diffusion model and corroborate phase statistics theory. 

PACS numbers: 42.25.Dd, 42.30.Ms, 61.43.Gt 
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I. INTRODUCTION 

Diffusion is one of the most widely encountered phe- 
nomena in physics. The dissolving of sugar in water, the 
transfer of heat in a wire and the transport of carriers in a 
photodiode are all examples of diffusion. These processes 
are all described by the same diffusion equation. This 
equation also describes the diffuse transport of waves in 
disordered scattering materials. An example of a diffus- 
ing wave is the transport of light through a cloud or a 
colloid suspension. Wave diffusion is not limited to light; 
acoustic waves, microwaves, quantum particles or even 
seismic waves behave completely analogously. 

The last couple of decennia wave diffusion has been of 
strong interest both from applied as well as fundamen- 
tal points of view. In contrast to classical particle dif- 
fusion, wave diffusion is influenced by interference. The 
recognition that phase plays an important role in wave 
diffusion forms the basis for applications like diffusing 
wave spectroscopy 1] and optical coherence tomography 
, which are invaluable tools in the analysis of colloidal 
systems and in the optical imaging of biological tissue. 
Fundamental interest is motivated especially by the par- 
allels between light diffusion and transport of electrons in 
mesoscopic systems. These parallels have been demon- 
strated by the observation of the optical equivalents of 
universal conductance fluctuations [3j and weak localiza- 
tion HQ. 

Multiply scattering media are characterized by the 
transport mean free path £ (the average distance a wave 
travels through the medium before becoming diffuse) and 
the diffusion constant D (the rate at which diffuse waves 
spread over the medium). For electrons t can be consider- 
ably smaller than the wavelength A of the electron. When 

^ ^/27r), electrons become localized and the diffusion 
constant vanishes 0, Q • This breakdown of diffusion is 
called Anderson localization. Anderson localization of 
microwaves has been observed in quasi- ID systems [8|. 
Observations at optical wavelengths Q, however, remain 
under debate 

Ell- 
in quasi-lD microwave experiments, localization was 
shown to have a distinct effect on the statistical distri- 
butions of the intensity § and the phase Recently 



it has become possible to perform dynamic electric field 
measurements also in the optical regime, which allow a 
study of the optical phase [l2l Il3| . These types of mea- 
surements provide a direct measurement of the phase of 
diffusing waves and they can give unambiguous proof of 
the presence of Anderson localization of light. 

Here we report our optical experiments that thor- 
oughly test wave-diffusion theory by measuring the am- 
plitude and the phase of light transmitted through a 
strongly scattering, non-localizing medium. Using the 
technique of ultrashort pulse interferometry [IJ, we 
have access to the time-resolved intensity, the frequency- 
resolved intensity and the statistical distribution of the 
phase delay time. We demonstrate five different ways 
of extracting the diffusion constant from this multitude 
of experimental data. By comparing the results of these 
five different methods, we test the diffusion model thor- 
oughly and moreover show how to interpret time-resolved 
and frequency-resolved measurements consistently. 

In Section [H] of this paper we present a model for dif- 
fusion through a slab. From this model we will derive 
both the frequency-dependent and the time-dependent 
behavior and identify characteristic parameters that can 
be extracted from experimental data. The setup for mea- 
suring both amplitude and phase of transmitted light is 
described in Section ITTT1 In Section [Tvl we present our 
results and devote special attention to the comparison 
of different techniques to measure the diffusion constant. 
Our conclusions are given in Section Ivl 



II. THEORY 

A. An exact solution to the diffusion equation 

We consider the diffusion of scalar waves through a slab 
of random material. The slab fills the space < z < L 
and is infinite in the other directions. In this geometry 
it is convenient to use Fourier transformed coordinates 
q± = (q x ,q y ) for the transverse directions. The slab 
is illuminated from the left (z < 0) by a pulse at time 
t = 0. Since the incident light quickly loses its direc- 
tionality due to scattering, it is possible to model the 
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incoming light by a diffuse source inside the material. In 
this paper we assume isotropic scattering. The inclusion 
of anisotropic scattering in the source function and the 
description of anisotropic diffusion are tremendous, and 
basically unsolved complications. At this point we use 
a source located at a depth zq ss £ 14] . Later we will 
use a more sophisticated source. Under these conditions, 
the ensemble averaged energy density of diffuse light / is 
described by the diffusion equation |15| . 

[d t - DV 2 + Da 2 ] I(q x ,z; t) = 6(z-z )S(t)S(q±). (1) 



In this equation a = y/3/(££ a ) is the absorption coeffi- 
cient corresponding to an absorption mean free path £ a . 
The right hand side of Eq. Q is the source term, where 
S , (qj_) describes the transverse distribution of the source 
and has the unit of energy. The total energy in the source 
pulse is given by S(q± = 0). 

The propagation of light is affected by the boundaries 
of the slab. It has been shown [0|ll3 that reflections at 
the surfaces impose mixed boundary conditions on the 
diffusion equation, 



0,J(qx,O;t) 
-d z I(q±,L;t) 



l(q±,0;t)/z el , 
I(q±,L;t)/z e2 , 



(2a) 
(2b) 



where z e i and z e2 are so called extrapolation lengths. 
In the diffusion model, their values are given by z e i,2 = 
2£{\ + i? 1)2 )/3(l - i?i, 2 ). The reflection coefficients i?i 
and i?2, correspond to the left and the right boundaries 
respectively. These coefficients can be estimated from 
Fresnel's law using the refractive indices of the dielectrics 
outside of the slab and the effective index of the random 
medium [l7| . 

We solve the diffusion equation (Eq. JIJ) with mixed 
boundary conditions analytically in the frequency do- 
main. This solution can conveniently be used to find 
the field correlation function, the total transmission and 
the average diffuse traversal time. We use the same ap- 
proach as in ^tJj with the exception that we extend the 
model to allow for different extrapolation lengths at the 
two boundaries and use an exponential distribution of 
the source intensity. 

When Eq. (JTJ is Laplace transformed with respect to 
t, an ex pres sion for the energy density / can be found 
directly 
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where we have defined r\ = ^JiQ/D + q 2 ± + a* 
Laplace transform parameter £1 describes the frequency 
of intensity oscillations and is much smaller than the op- 
tical frequency of the field u>. A and B are found from the 
boundary conditions 11' a I) and (|2b|) after tedious algebra 

«+(zo) - 2[z el n + l]e**° 



(3) 
The 



A( V ) 



B(r,) = 



7 



j-(L)e 



7 +(£~z )-2[z e2 ??+l]e"( L - z °> 
7-(i)e ? ' L 



(4) 
(5) 



where we defined function as 

7 ± ( a; ) = (z el r)+l)(z e2 v+l)e r > x ±(z el r)-l)(z e2 r 1 -l)e- r i x . 

We now have the exact solution to the diffusion equation 
with mixed boundary conditions. In order to find the 
transmitted intensity flux, we calculate the forward flux 
J z = —Dd z I at the slab surface z = L, 



7-(L) 



where F Zo is given by 



(7) 



(8) 



Eq. describes the transmission for a source located 
at depth Zq. A more realistic and more sophisticated 
model assumes an exponential distribution of the source 
light. The exponential distribution models how light be- 
comes diffuse by being scattered out of the incoming co- 
herent beam. We adapt Eq. Q for the exponential source 
model by convolving F Zo with a (normalized) exponential 
source function, 



r L cxp(-yi) 

W = / dZo £[l-eM-L/£)] FzMZQ) 
1 — cxp(Lrj — L/£) 1 + z e iTj 
~ 1 - exp(-L/f) l-£r] ' 



(9) 



We obtain the transmitted flux for the exponential source 
from Eq. JJJ) by simply replacing F Zo by Ft. 

An important quantity in the analysis of random me- 
dia is the total transmission. The total transmission is 
found by integrating the flux over the whole back surface 
of the sample (this corresponds to taking = 0) and 
integrating over time (fi = 0). The ensemble averaged 
total transmission coefficient Ttot is therefore defined as 



J z {q ± = 0,n = 0) J z {v = a) 



S(qx = 0) 



S(qx = 0) 



(10) 



Neglecting absorption, the equation is evaluated to re- 
produce the well known result |19| 



+ Zel 



L + z el + z e2 



+ O(exp(-L/0). (11) 



This relation between Ttot, £ and L is often used to de- 
termine the mean free path experimentally by varying 
L. 

Next, we calculate the electric field correlation function 
Ce for the transmitted light. This correlation function 
contains information about the dynamics of the diffusion 
process, 



C E (Q) 



(E(uj)E*(uj + ri)) _ J z (^iQ./D + a 2 ,L) 



(\E(u;)\)(\E(uj + m 



J z (a,L) 



(12) 
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where E(lo) is the complex field amplitude of the trans- 
mitted light for a incoming field of optical frequency w 
and unit amplitude. The brackets ( ) are used to ex- 
plicitly denote ensemble averaging over all possible con- 
figurations of the disordered sample. We obtained the 
right hand side by assuming ergodicity and applying the 
Wiener-Khinchin theorem. 

The field-field correlation function in Eq. 1|12[1 is the ex- 
act result for diffuse transport through a slab using mixed 
boundary conditions and an exponential source distribu- 
tion. Earlier results (Refs. 01 an d H3) are reproduced 
by using the simpler sheet source representation of Eq. 
(JSJ and choosing z e \ = z e 2 or z e \ = z e i = 0. 

We will now turn to the intensity correlation function. 
This function relates two single channel transmission co- 
efficients. The single channel transmission coefficient T 
describes transmission from one input angle to one out- 
put angle. Integrating T over all outgoing angles yields 
the total transmission coefficient T tot . The intensity cor- 
relation function is defined as 



Ci(Q) 



_ (5T(uj)5T(oj + fl)) 
= (TH> (T(u + fi)) ' 



(13) 



where ST(lo) = T (lu) — (T (lu)) . A well known approxi- 
mation for the intensity correlation function is given by 



c/(n) = \c E m' 



(14) 



Eq. 1I4J1 is referred to as the C\ approximation and is 
valid for diffusive transport in multiply scattering media 
far away from the localization transition |2l|. Eqs. i|12|) 
and H14H show that both Ce{Q) and C/(0) depend on 
the diffusion constant only by means of the reduced fre- 
quency Vt/D. Fitting the frequency dependence of C/ is a 
commonly used method to extract the diffusion constant 
from measured correlation functions. 

It is instructive to introduce the characteristic traversal 
time for diffusive transmission r t , which is defined as 
the average time it takes a pulse of light to travel through 
the medium, 



n 



J dt,J z (t,L)t 
J dtJjt, L) 



-i lim 



dc E {n) 
on ' 



(15) 



The right hand side was obtained by rewriting the defini- 
tion of Tt in the Laplace domain representation and using 
Eq. I|12|). For zero absorption we find 



L 2 - U 2 



Ti 



6D 



3L e D 



0{e- h l"). 

(16) 



The diffuse traversal time is of fundamental interest since 
it relates to the Thouless criterion for localization [23| . 
Furthermore, the time scale is of practical interest since 
measuring r t provides a method of determining the dif- 
fusion constant. Our result in Eq. (| 1 fc>|) gives corrections 
of order z e L/D and higher to the value of Tt = L 2 /6D 



found by Landauer et al. |22(. These corrections are es- 
pecially relevant when L/z e < 10, which is the case for 
thin samples or samples with a high extrapolation length 
due to internal reflection. 



B. Phase statistics 

The crucial difference between diffusion of particles 
and wave diffusion is interference. For this reason we are 
interested in the phase of light that propagates through 
a scattering medium. Analysis of phase information is 
complementary to the analysis of the intensity and pro- 
vides an independent method of measuring the traversal 
time t 4 and therefore the diffusion constant. We consider 
only single channel phase statistics, which means that 
we relate phase and amplitude for one input angle to the 
phase and amplitude for a single output angle. 

Since the diffusion equation only describes the average 
intensity, an extension is needed in order to predict phase 
statistics. The statistical properties of the phase were 
predicted by van Tiggelen et al. |24j by assuming Gaus- 
sian statistics of the transmitted field J33j • This Gaussian 
assumption is valid when a high number of independent 
paths contributes to the field at the back surface of the 
random material. The central limit theorem predicts that 
in this situation the real and imaginary parts of the fields 
are described by a normal distribution [25j . Equivalently, 
the field amplitude is Rayleigh distributed and the phase 
</> has a uniform distribution between and 2tt. Neither 
the distribution of the intensity nor the distribution of 
the phase contains information about the diffusion pro- 
cess. Much more interesting is the probability distribu- 
tion of the group velocity delay time 4>' = d<fi/duj. This 
probability distribution reflects dynamic properties of the 
diffusion process and provides a method of measuring 
the diffusion constant. The statistics of the delay time 
</>' were calculated in Ref. |24|. For this calculation the 
Gaussian field statistics were extended to describe the 
correlations of two fields at nearby frequencies. These 
correlations are given by the field-field correlation func- 
tion. The resulting joint Gaussian distribution was sub- 
sequently used to calculate the probability distribution 
of the delay time, 



Q 



i) 2 + Q 



3/2 ' 



(17) 



where cf>' = (/)'/((/)') and Q is a dimensionless parameter. 
(<j)') and Q can be calculated from the first and second 
terms in the Taylor expansion of the field-field correlation 
function: Ce — 1 + irt^l — Ml 2 + 0(f2 3 ), which results in 
{(j}') = n and Q = 26/ t 2 -10- 

In Ref. [3| the correlation function for a system with 
simplified boundary conditions was used to calculate r t 
and Q. Here it was shown that without absorption Q 
equals 2/5 while with absorption Q is reduced. However, 
by carefully examining our solution for mixed boundary 
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conditions, Eq. (|12(l . we find that Q increases above 2/5 
when the extrapolation lengths are nonzero. 

The intensity-weighted delay time W is a fundamental 
quantity since the sum of this quantity over all incoming 
and outgoing an gles equals 7r times the density of states in 
the medium [24l2fi| . The weighted delay time is defined 



W ee T<t>'. 



(18) 



T and 4>' are statistically dependent variables; for chan- 
nels with a low transmission the probability distribution 
of (j)' is broader [27]]. Because of the statistical depen- 
dency, the statistics of the weighted delay time cannot 
be deduced from the individual probability distributions 
of T and <fi' and has to be calculated on its own. The 
probability distribution of W was calculated in a similar 
way as the distribution of (j)' and is given by [24| : 



P(W) 



1 



VT+Q 



exp 



-2\W\ 



sgn(IU) + VTTQ 



(19) 



where sgn is the signum function[34| and W = W/ (W). 
The average weighted delay time was found [24| to relate 
to the diffuse traversal time according to (W) — (T) r t . 

The correlation function of the weighted delay time 
Cw is defined as: 



Cw 



(w(u)w(oj + n)) 

(W(uj)) (W(lo + Q)) 



(20) 



This correlation function was calculated in the C\ ap- 
proximation (Eq. CUP using a joint Gaussian distribu- 
tion that relates the fields at four frequencies 



C W (Q) 



1 

2^ 



dC E {tt) 



on 



Re[C E (nf ° B 
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Microwave experiments showed that deviations from the 
C\ approximation [C-x and C3 correlations) cause Cw to 
decay with frequency much slower than is described by 
Eq. I|21|) 27]. Therefore measuring Cw provides a good 
way of testing the validity of the C\ approximation and 
of looking for signs of localization. 



C. Diffusion in the time domain 

Although we found an exact solution to the diffusion 
equation in the frequency domain, the time-domain be- 
havior is not obvious from Eq. Q. In this section we 
analyze diffusion in the time domain and we present an 
alternative technique for finding the diffusion constant. 
The time-resolved transmission can, in principle, be cal- 
culated by inverse Laplace transforming Eq. J7|) by means 
of contour integration [2^] . Unfortunately Eq. (7J) has an 
infinite number of poles, none of which can be found an- 
alytically when the extrapolation lengths are non-zero. 



Using a different approach we will show that the diffu- 
sion constant can be found by analyzing only the long- 
time behavior of the transmitted flux. 

A complete set of solutions to the diffusion equation 
(Eq. Q) is given by: 

I qr .g{q ±1 z;t) = sm(q z z + 9) exp (-[q 2 + a 2 ]Dt) 9(f) 

where q 2 = q\ + q 2 and Q(t) is the Heaviside step func- 
tion. In an infinite medium the longitudinal spatial fre- 
quency q z and phase 9 can be chosen freely. In a finite 
slab, however, there is an infinite, discrete set of combi- 
nations of q z and 9 for which the boundary conditions are 
fulfilled. For the boundary conditions given by Eqs. I|2a|) 
and (|2 b|) . permitted values of q z and the corresponding 
9 can be calculated numerically. Every solution for q z 
corresponds to two poles in Eq. Q with r\ = ±iq z . 

We will only calculate the long-time behavior of diffu- 
sion. In the long-time limit only the solution with the 
lowest q z survives, since according to Eq. (|22|l all other 
solutions decay faster. We number this particular solu- 
tion q z i, 9\. Now we are able to calculate the diffuse flux 
for t » l/q 2 zl D, 

J z {q ± ,z;t) = -J (q±)cos(q z iz+6'i)exp(-[qf 1 + a 2 ]Dt) 

(23) 

Jo can be calculated by contour integrating Eq. Q) 
around the poles at 77 = ±iq z \. In this article we are in- 
terested only in the exponential decay time of the trans- 
mitted flux and therefore will not explicitly specify Jo- 
For the total flux (qj_ = 0) we find an exponential decay 
with a decay time Td, 



D, 



D, 



(24a) 
(24b) 



where the approximate solution in Eq. (|24b|) was found 
by linearly extrapolating 7(qj_,z;t) at the slab bound- 
aries (this is equivalent to the method of mirror images 
used in Ref. |2(Al29ip and L e = L + z e i + z e 2 is an effective 
slab thickness. The approximate solution (Eq. 124b!) ') can 
be used for thick samples (L 3> z e i, Z^a). 

It is interesting to notice the differences between the 
decay time and the diffuse traversal time Tt- The decay 
time Td describes the long-time decay rate of the energy 
density of diffuse light in the sample. This decay rate is 
given by the slowest term in Eq. 11' I'll and does not de- 
pend on the distribution of the source intensity. The dif- 
fuse traversal time, on the other hand, has contributions 
from all terms in Eq. (1221) and is mainly determined by 
the short-time transmission. The diffuse traversal time 
does depend on the distribution of the source intensity. 
Concluding, Td and r t are time scales that correspond to 
different aspects of diffusion. Therefore, the consistency 
of the diffusion model can be tested experimentally by 
measuring both Td and Tt for a series of samples. 
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D. Apparent non-exponential decay in a realistic 

experimental configuration 

In the previous section we found that the total trans- 
mitted flux decays single exponentially in the long-time 
limit. In an actual experimental geometry, however, it 
is not possible to collect all the transmitted light; only 
a finite area at the back surface of the sample can be 
imaged on the detector. We model the limited area by 
means of a Gaussian detection efficiency with a known 
waist Wd- Furthermore, we assume that the source light, 
<S(q_i_), has a Gaussian intensity distribution with waist 
w s . The total intensity reaching the detector Jd e t is found 
by integrating over all spatial frequencies q^, 

TTT/) f [ 

JdctW = J dq x J z (<i L ;L,t)exp{--c^wl), (25) 

The intensity profile at the sample surface is time de- 
pendent according to Eq. (1221) since modulations with a 
high spatial frequency q± decay faster than those with 
a low spatial frequency. Since we detect only the flux 
from a finite area, our detection efficiency is time depen- 
dent as well. We define r/ = (w^ + Wj)/8D, being the 
characteristic time scale for the time-dependent detection 
efficiency, and find the total detected flux from Eq. (|25p. 

JdctW = w y^ D) M<U. = 0;L,t). (26) 

t+Tf 

This equation shows that a finite detection area imposes 
a non-exponential envelope on the detected transmission 
and increases the detected decay rate. For thicker sam- 
ples the additional decay will be more pronounced since 
the diffuse decay, as described by Td, is slower. As a 
result, the diffusion constant found from a linear fit of 
In J(t) is structurally overestimated. Usually the pre- 
factor in Eq. (|26[) is omitted, corresponding to the as- 
sumption that the detection system collects light from 
a large area (wd S> w s ). The consequences of omitting 
this correction can, however, be significant: in our ex- 
perimental configuration the correction results in up to 
a 25% modification of the measured diffusion constant. 

Naturally, the finite-area correction given by Ea. l2tj|) 
equally applies for diffusion in the frequency domain. 
Unfortunately, it is inconvenient to apply the correc- 
tion in the frequency domain analytically. Therefore 
we use a numerical fast Fourier transform to correct the 
frequency-resolved transmission, Eq. Q , and all derived 
quantities. 

E. Five ways of measuring the diffusion constant 

In Sections III AIII Cl we presented methods to calcu- 
late the frequency correlations, the phase statistics and 
the transient behavior of light diffusing through a slab 
of randomly scattering material. Two important time 



scales were identified: the diffuse traversal time Tt and 
the exponential decay time Td- 

In our experiments we will test the consistency of the 
diffusion model quantitatively by extracting the diffusion 
constant from experimental data using five different tech- 
niques. If the model is valid, we expect all techniques to 
yield the same diffusion constant. This will, however, 
only be the case when the boundaries and the source 
intensity distribution are accounted for correctly. There- 
fore a comparison of the diffusion constants, measured 
using different methods, provides an excellent way of test- 
ing our diffusion model. 

Method I In the first method, the diffusion constant is 
found from the diffuse traversal time r t . The dif- 
fuse traversal time is obtained from time-resolved 
transmission using the definition in Eq. I|I5[) . Af- 
ter applying the finite-area correction, the diffu- 
sion constant is found by means of Eq. I|ltj[l ■ Since 
the transmitted intensity decays exponentially the 
value of Tt depends mainly on the transmission at 
short time scales. 

Method II The second method is to measure the decay 
time Td by fitting the long-time decay of the trans- 
mitted flux. Subsequently, Eq. (|24a|l is used to find 
the diffusion constant. Since Method II relies on 
the time-resolved transmission at long time scales, 
the parts of the data used in Method I and Method 
II are nearly independent. 

Method III In the third method, the intensity correla- 
tion function is extracted from frequency-resolved 
measurements. Fitting Eq. (|I4|I to the measured 
correlation function yields the diffusion constant. 

Method IV The fourth method relies on the measured 
optical phase and makes use of the statistics derived 
for the phase of diffuse light. When the field obeys 
Gaussian statistics (</>') equals the diffuse traver- 
sal time T t . Consequently, Eq. i|I 6|) can be used to 
extract the diffusion constant from the measured 
phase. Since Method IV only uses phase informa- 
tion and Methods III only uses the measured inten- 
sity, these two methods are fully independent. 

Method V In the last method the diffusion constant is 
extracted from measurements of the weighted de- 
lay time. With (W) / (T) = r t we find the diffuse 
traversal time. As in Methods I and IV we calcu- 
late the diffusion constant using Eq. (|I6|) . It has 
been shown |3(j that Method V is mathematically 
equivalent to Method I. Therefore, we will only use 
Method V to verify the consistency of our data pro- 
cessing. 

All together we now have five different methods of mea- 
suring the diffusion constant. A comparison of the results 
of these methods provides a thorough test of the diffusion 
model and the phase statistics. Furthermore it enables 
an unambiguous determination of the diffusion constant. 
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FIG. 1: Schematic representation of the Mach-Zehnder inter- 
ferometer used in the setup. The first beam splitter (BSl) 
divides the light between a reference arm and a signal arm. 
The light in the signal arm is focused on the sample by lens 
Ll. The transmitted speckle pattern is collimated by lens 
L2 and recombined with the reference beam at beam splitter 
BS2. Since the reference arm is a few millimeters shorter than 
the sample arm, the signal pulse does not overlap the refer- 
ence pulse temporally. Finally, aperture Dl selects an area 
that is smaller than the typical speckle size and polarizer PI 
blocks light with a polarization perpendicular to that of the 
reference beam in order to increase the signal to noise ratio. 
The beam containing the signal pulse and the reference pulse 
is propagated into a scanning interferometer (FTIR). When 
a sample is placed in the signal arm, only a fraction of the 
incident light reaches the FTIR. In order to balance the inter- 
ferometer, we use beam splitters (BSl and BS2) that reflect 
approximately 4%. 



III. EXPERIMENT 

We have presented a theoretical framework connecting 
time-resolved measurements to phase statistics and fre- 
quency correlations. In order to test this framework, we 
need to measure both the amplitude and the phase of the 
multiple scattered light over a range of optical frequen- 
cies simultaneously. We perform these measurements us- 
ing the technique of femtosecond pulse interferometry as 
described in Ref. [12j. This technique involves two in- 
terferometers. The first interferometer is of the Mach- 
Zehnder type and is shown schematically in Figure ^ A 
beam splitter divides the incoming light between a signal 
arm and a reference arm. In the signal arm the light of 
the laser is focussed on the sample to a waist diameter of 
approximately 30 /im using a lens with a focal length of 
6 cm. In order to probe different random configurations 
of scatterers we illuminate different areas of the sample 
by translating the sample perpendicular to the incoming 
beam. For every sample position the transmitted light 
forms a different volume speckle pattern. The speckle is 
collimated using a second 6 cm lens and an area smaller 
than a typical speckle spot is selected from the pattern 
using an aperture with a diameter of 0.8 mm. At the 
second beam splitter the light transmitted through this 
aperture is combined with the reference pulse yielding a 



beam with two temporally separated pulses. 

The double-pulsed signal is directed into a Fourier 
transform infrared interferometer (FTIR). The FTIR 
(Biorad FTS-60A) is a Michelson interferometer and 
scans the delay time between two copies of the signal. 
A detector directly behind the FTIR obtains the field 
autocorrelation function of the pulse pair as a function 
of the extra path length in the scanning arm of the in- 
terferometer. Because of the temporal separation of the 
signal and reference pulses, it is possible to isolate the 
cross-correlate C(t) of the signal pulse with the refer- 
ence pulse. In the frequency domain the cross-correlate 
is given by 



C(w) = \S((j)\ 2 H a ((j)H r {Lj)E(w), 



(27) 



where S(uj) is the spectrum of the incoming pulse, H s (uj) 
and H r {uj) are the transfer functions of the signal and ref- 
erence arm of the Mach-Zehnder interferometer, respec- 
tively and E{uj) is the transfer function of the sample 
that we wish to extract. In order to find the transfer 
function of the sample, the cross-correlate is measured 
with and without the sample consecutively. Dividing the 
two functions yields the complex transfer function E(lo) 
containing both the phase and the amplitude of the trans- 
mitted light. Now the time-resolved field transmission 
E{t) can in principle be obtained by means of an inverse 
Fourier transform. In practice, however, the bandwidth 
of the transfer function is limited by the bandwidth of 
the source pulse. Outside this bandwidth the measured 
transfer function is dominated by noise, therefore addi- 
tional filtering is required before calculating E(t). In 
our case the pulses are generated by a Ti: sapphire laser 
(Tsunami, Spectra Physics) operating at 775 nm with a 
bandwidth of about 6 nm. We use a Chebyshev filter for 
filtering in order to have a minimum effect of side lobes 
and a maximum time resolution. 

In the experiment the signal and the reference beams 
have to overlap both at the aperture and at the detector 
in order to cause an interference signal. This condition 
implies that both the direction and the position of the 
signal beam are fixed and, as a result, the detection is 
limited to light emitted from a small area of the sample 
surface. Based on the geometry of the setup we approxi- 
mate the detection area by a Gaussian curve with a waist 
of Wd — 10 /im. 

We perform the measurements on samples consisting of 
a layer of rutile TiC>2 particles with a diameter between 
150 nm and 290 nm that are deposited on a substrate of 
fused silica. The titania grains have a refractive index 
of approximately 2.8. The extrapolation lengths can be 
calculated from the effective refractive index n e ff of the 
medium ^tJ- The effective index can adequately be es- 
timated from Mie theory [3l]]. For our samples we find 
n e g =1.34 and the corresponding extrapolation lengths 
are z e \/t = 0.69 for the left boundary and z e iji = 1.71 
for the right boundary. We measured the total trans- 
mission as a function of sample thickness and found a 
transport mean free path of I = 0.97 ± 0.10 ^m by fitting 
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FIG. 2: Time-resolved intensity transmission for a 10.1 /im 
thick sample consisting of Ti02 grains. The observed non- 
exponential decay (solid line) agrees with the finite-area cor- 
rection, Eq. 1261 . over four decades. Theoretical curves are 
obtained from Eq. (J7J and the shape of the filter. The dot- 
ted line was corrected for the time-dependent detection effi- 
ciency due to focussing, using Eq. 1261 with w a = 15 /im and 
Wd = 10 /im as estimated from the experimental configuration 
and a fitted diffusion constant of D — 27.0m 2 s _1 . The dashed 
fine is the uncorrected curve for the same diffusion constant. 
The inset shows an example of the interference signal at the 
detector as a function of the delay length in the scanning in- 
terferometer. The average intensity transmission is obtained 
from 50 such measurements performed on different areas of 
the sample. 

the data to Eq. {HJ . 

Our samples range in thickness between 1.5 ± 0.3 /im 
and 18.0 ± 0.3 /im. Since the samples are on a substrate 
that is much thicker than the layer of titania, it is neces- 
sary to compensate for the extra delay in the substrate. 
In order to accurately determine the extra path length, 
we direct the light that is reflected from the substrate 
into the FTIR without repositioning the sample. The 
thickness of the substrate is deduced from the time delay 
between the reflections from the front and the back of 
the substrate. 

With the setup described in this section we are able to 
measure the complex transfer function of random media. 
Below we analyze these transfer functions in the time do- 
main, the frequency domain and by looking at the phase 
statistics. 



IV. RESULTS 

A. Time domain measurements 

First we consider the decay in time of a transmitted 
pulse. The inset in Fig.[2]shows the raw data obtained by 
measuring the cross-correlate at a single position of the 
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FIG. 3: Diffuse traversal time rt (triangles pointing upward), 
and scaled decay time Tdiv 2 /6 (triangles pointing downward), 
as a function of sample thickness. Td is found by fitting the 
exponential decay of J(t) (Method II); the error bars indi- 
cate the values found from fitting the first part of the decay 
(lowest value) and the last part of the decay (highest value). 
The dashed line is the theoretical value of Td- The diffuse 
traversal time is found from the time-resolved intensity trans- 
mission obtained by numerically evaluating Eq. (1151 (Method 
I). Without taking into account detection efficiency, theory 
(dotted line) predicts that the two time scales converge for 
thick samples. When the theory is compensated for the effect 
of a finite detection area (solid line), this convergence is lost 
in agreement with the experimental data. 



sample. For different positions of the sample the trans- 
mitted pulse is distorted differently. We extract E(t) for 
every measurement as described in Section IIIII and aver- 
age the corresponding intensities over 50 sample positions 
to obtain the normalized transmission J(i) = (\E(t)\ 2 ). 
Fig.[3shows the time-resolved intensity transmission J(t) 
for a 10.1 /im thick sample. We find that in the long-time 
limit J(t) has a nearly exponential decay for more than 
four decades. The measurements are fitted with the the- 
oretical curve obtained from Eq. J7J) convolved with the 
frequency filter that was used in the processing of the 
raw data. We find a good fit for a diffusion constant 
of D = 27.0 m 2 s _1 taking into account the effect of the 
limited area of detection. For comparison, the theoreti- 
cal curve without correction for the detected area is also 
shown in Fig. |3 The corrected curve exhibits a signifi- 
cantly faster decay, especially for t < 2 ps. 

In order to analyze the decay of J(t) more quantita- 
tively, we extract the diffuse traversal time r t and the 
decay time, Td from the measured flux. In Section llV Dl 
the diffusion constant will be calculated from these two 
times scales using Methods I and II (Section III Efl re- 
spectively. The first time scale r t is obtained from the 
time-resolved transmission directly using Eq. (|15fl . The 
second time scale Td is extracted from an exponential fit 
of the intensity decay. We fit the data between t > Td and 
the point where the intensity is dropped below the noise. 
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It was found that the time constant associated with the 
first half of this range is always significantly higher than 
the time constant for the second half, as is predicted by 
Eq. i|2t)|) . In Fig. [21 the decay time Td is compared to the 
diffuse traversal time r t . It was shown in Section [H] that 
the differences between the two time scales are caused by 
surface effects and the limited detection area. We find a 
good agreement with theory for the decay time as well as 
for the traversal time, both using the same fitted average 
diffusion constant of D = 25.5 m 2 s _1 . 



B. Phase statistics 

An independent way of measuring the diffuse traversal 
time is by analyzing the phase information. For different 
positions of the sample we obtain the phase <fi from the 
complex transfer function E(lj) that was measured using 
the technique described in Section ITTll For all different 
frequencies in the 6 nm bandwidth of the measurements 
we calculate the delay time <p'{w) = d(j)(uj) / duj and the 
weighted delay time W(uj) = cf>' (uj)\E(uj)\ 2 . By binning 
the values of </>' and W, we obtain the probability distri- 
butions shown in Fig. 0] The distributions are in good 
agreement with the predicted functional forms from the- 
ory, Eq. (|f 7(1 and Eq. I)IH[I. This agreement is a clear ex- 
perimental proof that the transmitted light is described 
well by a circular complex Gaussian distribution. For a 
sample with a thickness of lO.l^m, the characteristic pa- 
rameter Q determining the width of the distribution is 
calculated to be Q = 0.44. The experimental data gives 
Q = 0.47, corresponding to a slightly lower maximum of 
P{4>')- The high value of Q indicates that there is no 
measurable effect of absorption (which would decrease 
Q). Moreover, Q is clearly larger than the value of 2/5 
predicted in Ref. [24[. This observation shows that even 
for thick samples [L 10^) the effect of reflections at the 
surfaces cannot be neglected. 

We obtain the diffuse traversal time using r t = (cj)') 
(Method IV) and r t = (W) (Method V) and compare 
these results to the value found from the time-resolved 
intensity measurements (Method I). Fig. shows as 
obtained by these three different methods. The values of 
((/>') coincide almost perfectly with Tt found from time- 
resolved analysis. This agreement again confirms the ex- 
cellent validity of the C\ approximation and the theory 
of phase statistics. 

As expected, the results from Methods I and V agree 
very well. Although these methods are equivalent in the- 
ory, the time-domain data, on which Method I is based, 
have been filtered (see Section HID , whereas Method V 
uses the unfiltered frequency-domain measurements di- 
rectly. Since the differences between the values obtained 
by Methods I and V are minute, we conclude that the 
determination of r t is insensitive to frequency domain 
filtering. In Section IIV Dl we will use r t to extract the 
diffusion constant for each sample. 




4>7<4>*> 




FIG. 4: Probability distributions for the delay time cj>' (top) 
and the weighted delay time W (bottom) as measured in a 
10.1 /im thick Ti02 sample. The dimensionless parameter 
Q characterizes the width of these distributions. We find 
Q — 0.47 from a fit of the theoretical curves given by Eqs. 
1171 and 119H (solid lines). The average values of 4>' an d W 
are used in Methods IV and V respectively to find the diffusion 
constant. 



C. Frequency domain measurements 

In Sections llV Al and IIV Bl we presented measurements 
of the traversal time Tt and Q, the characteristic param- 
eter for phase statistics. These two parameters are re- 
lated to the first and second order terms in the Tay- 
lor expansion of the field-field correlation function Ce 
around Q = 0. In this section we go a step further and 
investigate the full frequency correlation functions of the 
transmitted light. We investigate the field-field correla- 
tion function, the intensity correlation function and the 
correlation function of the weighted delay time consecu- 
tively. 

We first look at the field-field correlation function. 
This function is related to the timc-resolvcd intensity 
transmission by a Fourier transform and provides an al- 
ternative way of studying the propagation of diffuse in- 
tensity without having to worry about possible artefacts 
introduced by filtering. In analyzing the time-resolved 
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FIG. 5: Diffuse traversal time Tt measured using three dif- 
ferent techniques. Method I (squares) calculates Tt from the 
time-resolved intensity. Method IV (triangles) obtains Tt from 
the measured optical phase alone, whereas Method V (circles) 
uses the intensity-weighted phase information. The excellent 
agreement indicates that the transmitted field is described 
by Gaussian distribution. The solid line are the theoretical 
values for a diffusion constant of D — 25.5 m 2 s -1 . 




FIG. 6: Measured field correlation function Ce (real part: cir- 
cles, imaginary part: squares) for a Ti02 sample of thickness 
10.1 fj,m as a function of the frequency difference Q between 
two optical frequencies. The horizontal axis is scaled by r t _1 
as found by measuring {<fi'}. Excellent agreement with the- 
ory (solid lines) confirms the diffusion model with boundary 
corrections. 



transmission plotted in Fig. [21 we only extracted two pa- 
rameters, T t and Td- Whereas the measured time-resolved 
transmission curve showed some minor fluctuations com- 
pared to theory, the field-field correlation function is per- 
fectly smooth up to il = 15t* and we find that the the- 
oretical curve matches the experimental data very well, 
as is shown in Fig. 

In order to test the C\ approximation (Eq. I|14|)l di- 
rectly, we examine the intensity correlation function 




FIG. 7: Measured correlation functions for the intensity 
Ci(yi) (diamonds), and the weighted delay time Cw(fi) (cir- 
cles), in a 10.1 fj,m thick Ti02 sample. The solid lines are the 
Ci approximations for both correlation functions. Except for 
some spurious oscillations, agreement with theory is evident. 



Cj(fi). The presence of long range (C2) and infinite range 
(C3) correlations would show up by comparing the inten- 
sity correlation function to the C\ contribution. In our 
experiment, however, we find a good agreement to the 
C\ theory as is shown in Fig. [7| At f2r t sa 8 a slight 
deviation of unknown origin is found in the correlation 
function. Surprisingly this deviation was absent in the 
field- field correlation function Ce ■ We determine the dif- 
fusion constant by fitting the intensity correlation func- 
tion and find a diffusion constant of 27 ± 3 m 2 s _1 for a 
10.1 /im thick sample. 

Finally, we present the correlation function for the 
weighted delay time CW(f2) in Fig. Apart from the 
same deviations that were found in the intensity corre- 
lation function, the agreement with Eq. l|21f) is evident. 
This observation provides the first experimental confir- 
mation of the correlation function of the weighted delay 
time at optical frequencies. 



D. The diffusion constant 

Altogether we have presented five different methods of 
determining the diffusion constant experimentally. Meth- 
ods I and II use the measured time-resolved intensity 
transmission to find two time scales, Td and r t , from 
which the diffusion constant can be calculated. Subse- 
quently, we showed that r t can be obtained from phase 
statistics in two different ways by analyzing the delay 
time (Method IV) and the weighted delay time (Method 
V). Finally, we measured the diffusion constant by fitting 
the intensity correlation function (Method III). The re- 
sults of these five methods are summarized in Fig. [S] for 
nine samples of different thickness. We find that all dif- 
ferent methods yield the same diffusion constant, within 
the experimental accuracy, for a given sample. This ob- 
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FIG. 8: Measured diffusion constant for samples of different 
thicknesses obtained in five different ways. The top figure 
shows the diffusion constant obtained from the diffuse traver- 
sal time Tt. We measured the traversal time directly in the 
time domain (Method I, squares), by extracting the average 
of the delay time (</>'} (Method IV, circles), and the weighted 
delay time (W) (Method V, triangles). Error bars are only 
presented for {(/)'}, but the errors in the two other values are 
comparable. The diffusion constants in the middle figure are 
calculated from the decay time Td which is obtained by fitting 
the decay of the transmitted flux (Method II). The bottom 
plot displays the diffusion constants obtained from fitting the 
intensity correlation function (Method III). Except for the two 
thinnest samples, the consistency of the five different measure- 
ment methods is evident. 



servation is an experimental proof of the consistency of 
the diffusion model that was presented in Section [D] An 
average diffusion constant of D = 25.5 ±1.0 m 2 s _1 is 
found; individual samples with different thickness have 
slightly different values of the diffusion constant ranging 
from 19m 2 s _1 to 28m 2 s~ 1 . Since these variations in the 
diffusion constants are reproduced in for all methods, we 
conclude that the scatter is a result of the varying sample 
structure and that it is not the result of a measurement 
error. 

The error bars in Fig. [H] are derived from the uncer- 
tainty in the sample thickness and the uncertainty in de- 
termining the actual value (r t , Td, or the curve fit to Cj) 



that was used to calculate the diffusion constant. For 
thinner samples we find larger error bars since the uncer- 
tainty in the thickness is relatively large. For the three 
methods that are based on the traversal time r t we find 
that the uncertainty in D decreases with increasing sam- 
ple thickness. The determination of the decay time r<j 
on the other hand, becomes increasingly more inaccurate 
because of the non-exponential decay of the transmitted 
intensity. 

The diffusion constants obtained by fitting the decay 
of the transmitted intensity (Method II) appear to be 
slightly lower for thinner samples. This behavior is con- 
sistent with earlier observations . It should be noted, 
however, that the decay time in these samples is com- 
parable to the time-domain resolution. Measurements 
based on phase information (which are not limited by 
the time resolution) show no thickness dependence of the 
diffusion constant. 

For the thinnest sample we find different diffusion con- 
stants depending on which method we use. The differ- 
ence is most apparent when comparing Methods I and 
V to Method IV. Both the direct time-domain measure- 
ment of Tt (Method I) and the average weighted delay 
time (Method V) are lower than expected, resulting in 
significantly higher values of the measured diffusion con- 
stant. The difference in observed diffusion constants is a 
clear indication that the transmitted field does not have 
a Gaussian distribution. The discrepancy between the 
different values of Tt is consistent with an increased trans- 
mission at short times. Therefore this observation sug- 
gests an influence of coherent transmission or single scat- 
tering. Especially, it shows that shorter traversal times 
are associated with higher intensities, since the low val- 
ues of Tt are not reproduced in the unweighted delay time 
measurements (Method IV). 



V. CONCLUSION 

Pulse interferometric measurements allow a sensitive 
determination of the complex transfer function of a ran- 
dom medium. This transfer function can be used to cal- 
culate the time-resolved intensity, the frequency-resolved 
intensity and notably the phase of diffuse light. These 
three complementary sets of data are obtained in a sin- 
gle measurement. 

In Section [H] we presented a consistent theoretical 
framework to interpret the experimental data. The 
framework is built around an exact solution to the dif- 
fusion equation with mixed boundary conditions. Our 
solution is a generalized version of the result presented 
in Ref. an d uses a more accurate description of the 
source intensity distribution. A new result is the non- 
exponential envelope that is imposed on the detected flux 
due to a finite detection area. This envelope contributes 
to the detected intensity decay and should be taken into 
account when extracting the diffusion constant. 

By analyzing the diffusion theory both in the time do- 
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main and in the frequency domain we have identified two 
relevant time scales. The first is the diffuse traversal time 
r t . The traversal time is the natural time scale in the 
analysis of phase statistics and correlation functions. The 
second time scale is the decay time Td associated with the 
exponential decay of diffuse intensity. These two times 
are affected differently by the boundary conditions and 
the finite area of detection. 

We measured the diffusion constant using five differ- 
ent methods, which allowed a comparison of time do- 
main measurements, frequency domain measurements 
and phase measurements. For our samples consisting of 
TiC>2 particles we found that the results of the five com- 
plementary techniques agree almost perfectly for samples 
thicker than twice the mean free path. Our observations 
are a strong experimental proof of the validity of the dif- 
fusion model and the phase statistics theory. 

Finally we measured the correlation function of the 
weighted delay time. To our knowledge, this is the first 
report of such a measurement at optical wavelengths. 
The experimental data agree very well with the measured 
intensity correlation function and the predictions from 
phase statistics theory. These measurements demon- 
strate that it is possible to record phase-related corre- 



lation functions at optical wavelengths. 

Our analysis clearly shows that care has to be taken 
in including proper boundary conditions and correct- 
ing for the detection efficiency even for samples much 
thicker than the mean free path. Provided these effects 
are taken into account properly, the model used to de- 
scribe propagation of light through a random medium 
is consistent for all different methods of analysis. Our 
experiment demonstrates the versatility and reliability 
of pulse-interferometric measurements and validates the 
use of phase-sensitive quantities for the identification of 
long range correlations and possibly localization of light. 
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